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Abstract 


Singularly perturbed robin type boundary value problems with discontinu- 
ous source terms applicable in geophysical fluid are considered. Due to the 
discontinuity, interior layers appear in the solution. To fit the interior and 
boundary layers, a fitted nonstandard numerical method is constructed. 
To treat the robin boundary condition, we use a finite difference formula. 
The stability and parameter uniform convergence of the proposed method 
is proved. To validate the applicability of the scheme, two model problems 
are considered for numerical experimentation and solved for different values 
of the perturbation parameter, c, and mesh size, h. The numerical result 
is tabulated, and it is observed that the present method is more accurate 
and uniformly convergent with order of convergence of O(h). 
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1 Introduction 


Singular perturbation problems model convection-diffusion processes in ap- 
plied mathematics that arise in diverse areas, including linearized Navier— 
Stokes equation at high Reynolds number and the drift-diffusion equation 
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of semiconductor device modeling, heat and mass transfer at high Pe’clet 
number, and so on; see [6, 13, 18, 19]. The novel aspect of the problem 
under consideration is that we take a source term in the differential equa- 
tion that has a jump discontinuity at one or more points in the interior 
of the domain. This gives rise to an interior layer in the exact solution 
of the problem, in addition to the boundary layer at the outflow boundary 
point. Problems with discontinuous data were treated theoretically, in the 
case of the solution of the convection—diffusion with Dirichlet case problem; 
see [9, 10]. Authors of [2, 8, 14] discussed a self-adjoint Dirichlet type problem 
with a discontinuous source term. Authors [14, 15, 17] have examined two 
parameter singularly perturbed boundary value problems for second-order 
ordinary differential equations with discontinuous source term. Authors of 
[5, 7] discussed fitted nonstandard finite difference methods for singularly 
perturbed second-order ordinary differential equations. Singularly perturbed 
delay differential equation was examined by Mohapatra and Natesan [12] on 
an adaptively generated grid. Recently, Shandru and shanthi [3] presented a 
fitted mesh method to solve singularly perturbed robin type boundary value 
problems with discontinuous source terms. Indeed, still, there is a room to 
increase the accuracy and show the parameter uniform convergence because 
the treatment of singular perturbation problem is not trivial distributions and 
the solution is pended on perturbation parameter, ¢ and mesh size, h; see [6]. 
Due to this, the numerical treatment of singularly perturbed boundary value 
problems is need improvement. Therefore, it is important to develop a more 
accurate and convergent numerical method for solving singularly perturbed 
boundary value problems under consideration. 


2 Definition of the problem 


Consider the following singularly perturbed problem with Robin boundary 
condition of the form 


Ly(a) = ey” (x) + a(x)y'(x) — W(a)y(@) = f(a), ve Q7UQt. (1) 
Subject to boundary conditions 


oun = a1y(0) — Brey’( 
Lyy(1) = agy(1) + Boy’ 


where a1,0; > 0, a1 + 61 > 0, a2 > 0,82 > 0, and « > 0 is a small 
parameter. The functions a(x) and b(x) are smooth on Q, such that a(x) > 
a > 0 and W(x) > b > 0. Furthermore, the notations for the domain are 
Q = (0,1), Q- = (0,d), and N* = (d,1), where d € Q stands for the jump 
in the source function. Boundary value problem of the governing problem 
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under consideration is a model confinement of a plasma column by reaction 
pressure and geophysical fluid dynamics; see [4]. 

The solution y(z) of (1)—-(2) has a boundary layer near x = 0 due to the 
perturbation parameter, ¢ and interior layer due to the discontinuous source 
term. 


3 Properties of continuous solution 


The differential operator for (1) is given by 


and it satisfies the following minimum principle for boundary value problems. 
The following lemmas [6] are necessary for the existence and uniqueness of 
the solution and for the problem to be well-posed. 


Lemma 1 (Continuous minimum principle). Suppose that the function 
y € C1(Q) NC? (Q-UA*) satisfies Liy(0) > 0, Ley(1) > 0, and Ly(x) < 
0, for all e €Q7 UNF and [y’](d) < 0. Then, y(x) >0 for alle eQ. 


Proof. For the proof, we refer to [3]. 


Lemma 2 (Stability result). Consider the boundary value problem (1)-(2) 
subject to the conditions a(x) > a > 0 and b(x) > b> 0. Ifye C1(Q)N 
C?(Q-UNT), then 


ll vy lax Cmax{| Liy(0) |, | Ley(1) |, | Ly le-ua+}- 


Proof. For the proof, see[3]. 


Lemma 3. For each integer k satisfying 0 < k < 4, the solution of (1)—(2) 
satisfies the bounds || y“* lay gay S Ce, 


Proof. For the proof, see [3]. 


Lemma 4. Let y- be the solution of (P-). Then, for k = 0,1, 2,3, 


| (x) |< C1 +e7* exp(—2)) for all x € (0, []. 


Proof. For the proof, see [1]. 
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4 Formulation of the method 


The theoretical basis of the nonstandard discrete numerical method is based 
on the development of the exact finite difference method. The author of [11] 
presented techniques and rules for developing nonstandard finite difference 
methods for different problem types. In Mickens’s rules, to develop a dis- 
crete scheme, the denominator function for the discrete derivatives must be 
expressed in terms of more complicated functions of step sizes than those used 
in the standard procedures. These complicated functions constitute a general 
property of the schemes, which is useful while designing reliable schemes for 
such problems. 


For the problem of the form in (1)—(2), in order to construct the exact 
finite difference scheme, we follow the procedures used in [1]. 
Let us consider the following singularly perturbed differential equation of the 
form 


ey" (x) + a(x)y"(x) — b(x)y(x) = f(a). (3) 

The constant coefficient homogeneous problems corresponding to (3) are 
ey!" (x) + ay’ (a) — by(x) = 0, (4) 
ey" (x) + ay'(x) = 0, (5) 


where a(x) > a and b(x) > b. Two linear independent solutions of (4) are 
exp(A;x) and exp(A2x), where 


—at Va? + 4eb 
M12 = ——— ee (6) 
E 

We discretize the domain [0,1] using the uniform mesh length Ax = h such 
that ON = {x; = 29 + ih, 1,2,...,N, 29 =0,tn =1,h = 7}, where N de 
notes the number of mesh points. We denote the approximate solution to y(z) 
at the grid point 2; by Y;. Now our main objective is to calculate the differ- 
ence equation, which has the same general solution as the differential equation 
(4) has at the grid point x; given by Y; = Ai exp(\12;) + Az exp(A22;). Using 
the theory of difference equations and the procedures used in [1], we have 


Yj-1 exp(A12i-1) exp(A22i-1) 
det | Y; exp(Aia;)) exp(Ae2;) | =0. (7) 
Yi-1 exp(Ai2i41) exp(A2zi+1) 


Simplifying (7), we obtain 


h ha? + 4eb h 
52) Vin1 + 2cosh( = =” yY; — exp( 5a Yin = 0, (8) 


Expl 
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which is an exact difference scheme for (4). 
After doing the arithmetic manipulation and rearrangement on (8), for the 
constant coefficient problem (5), we get 


oe gait —¥Y; 
BE (exp(%) — 1) h 


=0. (9) 
h h 

The denominator function becomes W? = as (exp (=) - 1). Adopting 
a 


this denominator function for the variable coefficient problem, we write it as 


ier (exp (=) 1) , (10) 
ay 


where W? is the function of ¢, a;, and h. 
By using the denominator function W? in to the main scheme, we obtain the 
difference scheme as 


Vier — 2Y; + Yi-1 
w2 


ai uf biYi = fi. (11) 


Ny — 
LY, =€ 
This can be written as three term recurrence relation as 


BYi1+hiY4+ GiYinr =H, 1=1,2,...,N—-1, (12) 


where Eo=@ B= - Fh, Gi= get F, and H = fi. 
To treat the boundary condition, we use the forward finite difference formula 
for 2 = 0 and the backward difference formula for i = N, respectively, for the 
first derivative term. 
That is, for 7 = 0, from (2), we have ay y(0)—fiey4 = p implies ay yo—fiey) = 
p, which yields 

Bre 


E 
(a1 +E )yp — PE, = p, (13) 


Similarly, for i = N, from (2), we have agy(N) + Soy'y = ¢ implies a2yn + 
Boy'y = q, which yields 


Bo Bo 
Pe ijeg = tie ye 14 
(a2 + >-)un — [-yn-1 = 4 (14) 
Therefore, Equation (1) with the given boundary conditions (2), can be solved 
using the schemes in (12), (13), and (14) which gives the N x N system of 
algebraic equations. 
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5 Uniform convergence analysis 


In this section, we need to show that the discrete scheme in (12) satisfies 
the discrete minimum principle and uniform convergence. Let us define the 
forward, backward, and second-order central finite difference operators as 

Yj+i 7 Yj 


¥j-¥j-1 D+Y; - D7¥j, 


DtY; = ; 


y= oy eee ae 
Lemma 5 (Discrete Minimum principle). Let V; be any mush function that 
satisfies v9 > 0, vy > 0, and Lv; < 0,1 =1,2,...,N—1. Then v,; > 
0,7=1,2,...,N. 

Proof. The proof is obtained by contradiction. Let j be such that V; = 


min V;, and suppose that V; < 0. Clearly, 7 ¢ {0,N}, Vj41 — Vj > 0, and 
V; — V;-1 < 0. Therefore, 


E aj 
LV; = 53 (Vina — 2Vj + Via) + (Vita — Vj) — bY; 
E Qi 
= palin — Vi) — (Vi — Via) + F (Vina — Vi) — BV5 
J 
>0 


— 


where the strict inequality holds if Vj; —V; > 0. This is a contradiction and 
therefore V; > 0. Since 7 is arbitrary, we have V; > 0, 1=1,2,...,N. 


We proved above that the discrete operator L’ satisfies the minimum 
principle. Next, we analyze the uniform convergence analysis. 

Using the Taylor series expansion, the bound for y(a#;-1) and y(a;41) at 
x; are as 


We obtain the bound for 


eens. <Cly"(aa)|, (15) 
ly" (xi) — Dt D~y(ai)| < Ch? |y (ai)]. 
Similarly, for the first derivative term, we have 

ly’ (wi) — DT y(wi)| < Chly (ai)I, (16) 


where ly) (x,)| = SUDz, €(x0,2N) ly (xi)I, k = 2,3,4. 


Theorem 1. Let the coefficients functions a(a) and the source function f(z) 
in (1)-(2) of the domain Q be sufficiently smooth, so that y(x) € C*[0, 1]. 
Then, the discrete solution Y; satisfies 
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exp Tari 
(yi —¥)|<h(1+ sup (—) | 
x€(0,1) E 


Proof. We consider the truncation error discretization as 


IZ™ (ys — ¥i)| =|L% yi — LYY;|, 


D*D7~h? 
<Cley; + aii, — {e P ys + a;DTy;}I, 


7 


DD 


au 


yi) taily; — Dt y)I, 
h2 
<Cely’ — Dt D-y;| + Cel( Fa —1)D*D-y;| + Chly? |, 
<Ceh?|y |+Chlyf| + Chlyi'|, 
<Ceh|y\| + Chlyi'|. 
We use the estimate ele —1| < Ch, which can be derived from (10). 


Ah 
Indeed, define p= “.p € (0,00) .Then, 


1 
exp(p) — 1 


By simplifying and writing the above equation explicitly, we obtain 


he 1 
ela — 1] =a,h| - - =: a;hQ(p). 


_ exp(p) — p—1 
(= expla) — 1) 


and we obtain that the limit is bounded as 


lim Q(p) = 5+ lim Q(0) =0. 


p—0 2’ 


Hence, for all p € (0,00), we have Q(p) < C. So, the error estimate in the 
discretization is bounded as 


ILN (ys — ¥i)| < Ceh?|y| + Chly/'|. (17) 


From (17) and boundedness of derivatives of solution in Lemma 4, we obtain 


(1+2~texp (=2")) 
€ 
+er| (1 4ePex (=S*)) | 


|Z (y(as) — Yi)| <Ceh? 
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(:+2%en (=F) 
e+e “exp 
€ 
rer] (1 4ePem (=S*)) | 
€ 
<Ch{1+ = sup (—) ; 
x€(0,1) E 


Most of the time during analysis, one encounters with exponential terms 
involving divided by the power function in ¢, which are always the main 
cause of worry. For their careful consideration while proving the ¢-uniform 
convergence, we prove the following lemma. 


<Ch? 


since e€~? > e7?, 


Lemma 6. For a fixed mesh and for ¢ > 0, it holds 


(= a) =0 nas 
em ’ reek aed cei | 


lim max 
e301<i<N-1 


exp(—2U=2i)) 
lim max ——=——_ ]=0,  m=1,2,3,..., 
e301<i<N-1 em 


where x; = th,h = 7,i=1,2,...,.N—1. 


Proof. Consider the partition [0,1] := {0 = 2% < a1 <-+++<ay_-1<a@y = 
1}. For the interior grid points, we have 


1<i<N-1 em = em em 
—a(1— 2;) —a(1— @y_-1) —ah 
exp exp exp | —— 
E E 
max < = : 
1<i<N-1 em em em 


as 41 =1—ay_1= h. 


—ah 
— (=) rm m! 
: = 0. 


lim ——~—+ = li ——— = li ——— 
e—>0 em ere exp(ahr) = aan (ah)™ exp(ahr) 


Theorem 2. Under the hypothesis of boundedness of discrete solution (i.e., 
it satisfies the discrete minimum principle), Lemma 6, and Theorem 1, the 
discrete solution satisfies the following bound: 
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sup max|y; — Y;|< CN7". (18) 
O<e<1l ? 


Proof. Results from boundedness of solution, Lemma 6, and Theorem 1 give 
the required estimates. 


6 Numerical examples and results 


To validate the established theoretical results, we perform numerical experi- 
ments using the model problems of the form in (1)—(2). 


Example 1. Consider the following problem: 


oo y(2)= f(x), ren unt, 
y(0)—ey’(0)=1, yQ)-y'0) =-1, 
where 
0.7, O<2<0.5, 
w= {h O5<a<1. 


Example 2. Consider the following problem: 


ew + tet (2) =f(x), rEeQ UN, 
y(0) —ey(0)=1, yQ)-yG) =1, 


where 


lta, O0<2<0.5, 
f(x) = 
4, 0.5<a<1. 


Having y; = yy (the approximated solution is obtained via the fitted operator 
finite difference method) for different values of h and ¢, the maximum errors. 
Since the exact solution is not available, the maximum errors (denoted by 
EN) are evaluated using the double mesh principle [6], for fitted operator 
finite difference methods using the formula 


EN := max _ 
: 0<j ax |yp v2; 


Furthermore, we will tabulate the ¢-uniform error 


EN = max EN. 
O<e<l 


The numerical rate of convergence is computed using the formula [6] 


AN. log(B&) — log(B2) 
e log(2) 
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and the e-uniform rate of convergence is computed using 


n _ log(E™) — log(B?") 


a log(2) 


Table 1: Maximum absolute errors for different values of ¢ and number of mesh size, 
N for Example 1. 


é N=32 N=64 N=128 N=256 N=512 
10-* ~—2.0313e-02 1.0156 e-02  5.0781e-03 = 2.5391e-03 — 1.2695e-03 
10-8 —2.0313e-02 1.0156 e-02 5.0781e-03 2.5391e-03 — 1.2695e-03 
10-'* —_2.0313e-02 1.0156 e-02 5.0781e-03 2.5391e-03 —1.2695e-03 
10~1® 2.0313e-02 1.0156 e-02 5.0781e-03  2.5391e-03  1.2695e-03 
10~°  2.0313e-02 1.0156 e-02 5.0781e-03 2.5391e-03 —1.2695e-03 


EN — 2.0313e-02 1.0156 e-02  5.0781e-03  2.5391e-03  1.2695e-03 
RN 1.0001 1.0000 1.0000 1.0001 


—— Ns32 


Numerical Solution 
Numerical Solution 


0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 06 0.7 08 09 1 
x 


Figure 1: Behavior of numerical solution at ¢ = 2-5 and different values of N for 
Examples 1 and 2, respectively. 


Table 2: Comparison of maximum absolute errors and order of convergence for Example 
1 at number of mesh points N. 


E N=64 N=128 N=256 N=512 N=1024 
Present method 
EN 1.0156e-02 5.0781e-03 2.5382e-03 1.2467e-03 5.5910e-04 
RN 1.0000 1.0000 1.0257 1.0257 
Method in [3] 
EN 2.5658e-02 1.4128e-02 7.5359e-03 3.8225e-03 1.7536e-03 


RN 0.8605 0.90671 0.9793 1.1242 
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Table 3: Maximum absolute errors for different values of ¢ and number of mesh size, 
N for Example 2. 

E N=32 N=64 N=128 N=256 N=512 
10-4 8.3431e-02  4.1854e-02 2.0962e-02 1.0489e-02 5.2329e-03 
10-8 = 8.3431e-02  4.1854e-02 2.0962e-02 1.0489e-02 5.2329e-03 
10-7 -8.3431e-02 4.1854e-02 2.0962e-02 1.0489e-02 5.2329e-03 
10-'6 =8.3431e-02 4.1854e-02 2.0962e-02 1.0489e-02 5.2329e-03 
10-79 8.3431e-02 4.1854e-02 2.0962e-02 1.0489e-02 5.2329e-03 


En 8.3431e-02 4.1854e-02 2.0962c-02 1.0489e-02 5.2329e-03 
eal 0.9952 0.9976 0.9989 1.0032 


Table 4: Comparison of maximum absolute errors and order of convergence for Example 
2 at number of mesh points N. 


€ N=64 N=128 N=256 N=512 N=1024 
Present method 
EN 4.1854e-02 2.0960ce-02 1.0408¢-02 4.7906ce-03 2.0192e-03 
RN 0.9977 1.0099 1.1194 1.2464 
Method in [3] 
EN 9.6698e-01 5.8056e-01 3.2795e-01 1.7313e-01 8.1501e-02 
RN 0.7364 0.8239 0.0.9216 1.0870 


0.014 0.04 
—— N=32 —— —— N=32 
—+— N=64 6.096. —+— N=64 || 

0.012 —o— Ne128 1] —o— N=128 


0.1 0.2 0.3 0.4 0.5 0. 
x 


Figure 2: Pointwise absolute error plot at ¢ = 2~° and different values of N for 
Examples 1 and 2, respectively. 


7 Discussion and conclusion 


This study introduced a uniformly convergent numerical method based on 
nonstandard finite difference method for solving singularly perturbed second- 
order ordinary differential equations of Robin type boundary value problems 
with discontinuous source term. Due to discontinuity in the source term, 
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— o(1/N) 
+ 104 
e108 
a 10-1? 

5 s Oe a10'® 

5 5 10° 

8 8 

2 2 

£ 2 £ 2 

10 § 10 

8 & 

= = 

10 10° 


10° 10° 10! 10° 10° 
N- Number of mesh points N- Number of mesh points 


Figure 3: e-uniform convergence with NSFDM in Log-Log scale for Examples 1 and 2, 
respectively. 


there is an interior layer occurring. To fit the interior and boundary layer, a 
suitable nonstandard finite difference method on uniform mesh is constructed. 
The numerical results are tabulated in terms of maximum absolute errors, 
numerical rate of convergence, and uniform errors (see Tables 1-4) and com- 
pared with the results of the previously developed numerical methods existing 
in the literature (Tables 2 and 4). Furthermore, to see the position of the 
boundary layer, we plot the behavior of the numerical solution (see Figure 
1), as the number of mesh points increases, the maximum pointwise errors 
decrease (see Figure 2) and the e-uniform convergence of the method was 
shown using the log-log plot (see Figure 3). Unlike other fitted operator fi- 
nite difference methods constructed in standard ways, the method that we 
presented in this paper is fairly simple to construct. 
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